Complete realization of energy landscapes and non-equilibrium trapping dynamics in small spin glass and optimization problems

Energy landscapes are high-dimensional surfaces underlie all physical systems, which determine crucially the energetic and behavioral dependence of the systems on variable configurations, but are difficult to be analyzed due to their high-dimensional nature. Here we introduce an approach to reveal for the complete energy landscapes of spin glasses and Boolean satisfiability problems with a small system size, and unravels their non-equilibrium dynamics at an arbitrary temperature for an arbitrarily long time. Remarkably, our results show that it can be less likely for the system to attain ground states when temperature decreases, due to trapping in individual local minima, which ceases at a different time, leading to multiple abrupt jumps in the ground-state probability. For large systems, we introduce a variant approach to extract partially the energy landscapes and observe both semi-analytically and in simulations similar phenomena. This work introduces new methodology to unravel the energy landscapes and non-equilibrium dynamics of glassy systems, and provides us with a clear, complete and new physical picture on their long-time behaviors inaccessible by existing approaches.

Energy landscapes of physical systems are high-dimensional surfaces representing the dependence of system energy on variable configurations.Similarly, cost or fitness landscapes can be defined for optimization problems.Their characteristics determine crucially the emergent behavior of these systems.For instance, spin glasses and ferromagnetic spin systems are believed to have energy landscapes with and without a large number of local minima respectively 1,2 ; a similar analogy is made with the algorithmic-hard and -easy phases of combinatorial optimization problems 3 .A way to unravel and analyze the complete energy landscape is thus crucial to our understanding of these glassy systems.
Nevertheless, even for small systems, revealing completely their energy landscapes is difficult since they are high-dimensional functions.Existing approaches often omit some features of the landscapes for a feasible characterization.For instance, disconnectivity graphs (DG) connect attraction basins in the state space and show hierarchically how they are repeatedly segmented into smaller basins as energy decreases 4 .DGs have been applied to analyze energy landscapes of systems from protein folding to machine learning [4][5][6] , and can be improved using principal component analyses 7 .However, DGs only show the segmentation into basins, without showing their entropy nor how states are exactly connected, especially as basins may have multiple instead of one connection to other states.Another common approach is multi-dimensional scaling (MDS), which aims to preserve the high-dimensional distance between two states in a plot of reduced dimension 8 .For instance, one may preserve the distance between pairs of states in one-dimensional plots 9 .Nevertheless, MDS only shows pairwise distance while dynamics on these systems are definitely more complex than pairwise interactions.
On the other hand, efforts have been made to reveal the relations among the ground states but omitting other parts of the energy landscape, which is a goal different from the present study.For instance, the hierarchical structure among ground states was revealed by examining the similarity in spin domains [10][11][12] .For constraint satisfiability problems, 13 showed that solutions are grouped into separate clusters, while similarity among solutions are examined through relaxing discrete variables to be continuous during optimization 14 .Since this line of research primarily investigates ground states but not the rest of the energy landscape such as higher-lying local minima, and since the space spanned by the ground states is much smaller than the whole configurational space, they can investigate larger systems compared to the present study which aims to reveal the whole configurational space.
In this study, we introduce an approach to reveal for the complete energy landscape of complex disordered systems such as spin glasses and Boolean satisfiability problems.The approach is feasible on small systems, while for larger systems, we introduce a variant approach to obtain a partial energy landscape.Here, we do not aim to contribute to the extensive studies in the so-called thermodynamics limit, i.e. systems with an infinite size, of spin glasses and Boolean satisfiability problems; instead, our goal is to offer insights on the similarity and difference between the behaviors of small systems and those theoretical predictions in the thermodynamic limit, by visualizing and analyzing the complete energy landscape in small systems, which is much less explored in existing studies.
The obtained energy landscapes allow us to compute semi-analytically the approximate non-equilibrium dynamics at an arbitrary temperature for an arbitrary long time, out of reach by simulations limited by modern computational capability.In contrary to our common belief, we show that it can be less likely for the systems to attain ground states when temperature decreases, due to trapping in local minima; as time increases, trapping in individual minima ceases at different time, leading to multiple abrupt jumps in the ground-state probability.Our findings also provide insights on the effectiveness of simulated annealing compared with fixed-temperature dynamics, whereas only an extremely long annealing process that allows an escape from local minima may guarantee a ground state 15 .All in all, our approach opens up a new platform for analyzing the non-equilibrium dynamics of glassy systems, and provides us with a clear, complete and new physical picture on their long-time behaviors inaccessible by existing approaches and numerics.

Methods
We consider a system with N Boolean variables s i = ±1 , such that i = 1, . . ., N and s denotes the N-tuple (s 1 , s 2 , . . ., s N ) representing a variable configuration.We then denote the energy or objective function of the system as E( s) .Here, we examine two glassy systems as examples, namely (i) spin glasses 1 and (ii) K-satisfia- bility problems 16 .Here, we studied random instances; the difference between the energy landscapes of random instances and those instances with planted solutions 17,18 is worth studying and will be studied elsewhere.
For spin glasses, each s is a configuration of Ising spins and E( s) is given by where J ij = +1 with a probability f + and otherwise J ij = −1 ; the adjacency matrix a ij = 0, 1 characterizes different graph topologies.We multiply E by a factor of 1/2, such that a single spin flip leads to a unit change in energy.Depending on the topology, the parameter f + and the temperature, the spin system exhibits various phases such as paramagnetic, ferromagnetic and spin glass phases 1, 19 .
For K-satisfiability problems, or K-Sat for short, we introduce M clauses of the form (s µ 1 ∨ s µ 2 ∨ • • • ∨ s µ K ) labeled by µ = 1, . . ., M , each with K variables or their negation; the symbol ∨ corresponds to the "or" logical relation and the variables with an overline are negated.In this case, E( s) is given by where randomly drawn J µ,k = ±1 corresponds to the presence of the original or the negated k-th variable in clause µ .With the factor of 1 2 K , each violated clause increases the energy by 1 and the total energy is equivalent to the number of violated clauses.The ground state of the system is attained when E = 0 , i.e. all clauses are satisfied.Depending on the ratio α = M/N , the system exhibits various phases including a satisfiable phase at small α with an algorithmic-easy and -hard regime, followed by an unsatisfiable phase at large α 20 .The phase transition between the "easy" and the "hard" phases is only well defined for systems in the thermodynamics limit, but not for small systems studied here; nevertheless, we show cases with α = 1 and α = 4 which correspond to the "easy" and the "hard" phases in the thermodynamics limit.
Since there are N Boolean variables in the above systems, there are 2 N different variable configurations.If we consider two configurations s a and s b to be connected in the configurational space if their hamming distance is |� s a − � s b | = 1 , i.e. they differ only in the state of a single variable, the configurational space is effectively an N-dimensional hypercube.
To present this hypercube as an energy landscape, we take advantage of the integer disorders J defined in the above systems, which lead to discrete energy levels.We then represent each variable configuration s as a node in a network; two nodes are connected by a link if their hamming distance is 1.Next, we arrange nodes with the same energy on the same horizontal level in the network, with lower-energy configurations located at a lower row.We call this the full energy landscape (FEL).For the sake of clear illustration, we show an example of FEL of a small 3-Sat toy problem with N = 5 and α = 4 in Fig. 1a.One can see clearly how the 2 N = 32 different configurations are connected and arranged in different energy levels.Nevertheless, as N further increases, the number of states increases exponentially and FELs quickly become computationally infeasible and difficult to be clearly visualized.
To simplify the energy landscape, we group connected nodes on the same energy level into clusters; we denote C to be the total number of clusters.Two clusters a and b are connected if any pair of their constituent variable configurations are connected; the weight w ab of the connection is the total number of links between their constituent configurations.We call this energy landscape the coarse-grained energy landscape (CEL).The corresponding CEL of the FEL in Fig. 1a is shown in Fig. 1b, where the number of nodes is reduced from One may remark that two configurations differ by a single spin or variable flip are usually considered indistinguishable in the theoretical analyses under the thermodynamics limit, in such case only energy barriers of a height of at least the order of O( √ N) is relevant.Nevertheless, the objective of this study is not to analyze small systems as in the theoretical analyses in the thermodynamic limit, but instead we do the opposite to examine the trace of insights from the predictions in the thermodynamic limit to finite systems via visualizing the complete energy landscapes.In small systems, such a single spin or variable flip is important, for instance, in simulated annealing a single spin flip leading to higher or lower energy does affect its transition probability.Thus, in our case of small systems, we define configuration clusters using single variable flips.

CELs with local minima
We first analyze the CEL of small K-satisfiability problems with K = 3 and α = 4 , i.e. when the system is in the so-called "Hard-Sat" regime in the thermodynamics limit, usually characterized by complex energy landscapes as we understand.Here we can only study systems up to a small size N = 25 , which may seem small compared to the studies on ground states [10][11][12] , but indeed the number of configurations visualized in our study is 2 25 ≈ 3.3 × 10 7 , which is much larger than the number of ground states, i.e.O(10 2 ) or O(10 3 ) , analyzed in these studies.
As shown in Fig. 2a, as N further increases, the ratio of the number of clusters in CEL to the total number of variable configurations, i.e.C/2 N , decreases exponentially with N, implying that an extensive number of configurations can be grouped in CEL for a clear presentation.This also implies that C ∝ 2 γ N , with γ < 1 .As shown in Fig. 2b, γ increases with α , implying that the structure of energy landscape is more complicated at large α , consistent with our understanding of algorithmic-hard regimes in the thermodynamic limit, compared with easy ones at small α .One can also see that γ approaches ln 2 as α increases, implying that clusters are increas- ingly composed of individual nodes as there are more distinct energy levels, consistent with the segmentation of solution space in studies which focuses on analyzing ground states 3 .
Interestingly, as shown in Fig. S1a of the Supplementary Information (SI), the exponent γ for K-Sat problems with different K and α collapses onto a common function of α/K 2 .This implies that the decrease of nodes by grouping configurations in CEL is universal for different values of N, M and K, which is further shown by the ratio (C/2 N ) 1/(1−γ ) collapsed onto a common exponential decay against N in Fig. S1(b).
Other than a large reduction in the number of nodes, another advantage of CEL is the identification of local minima.Since connected configurations with the same energy are grouped in clusters in CELs, one can easily identify local minima as clusters where all neighbors are of higher energy; such identification is not trivial in FEL since it is difficult to examine if there exists a path from a configuration to a lower-energy one without passing through higher-energy configurations.In the CEL in Fig. 1b, one can see that there exist two local minima (triangles) with E = 1.
We show in Fig. 3a the low-energy portion of another examplar CEL of spin glasses on random regular graphs (RRG) with N = 15 and f + = 0.5 ; since the configurations s and −� s have the same energy according to Eq. (1), one can observe a symmetric structure in the energy landscape as expected, including a pair of local minima at E = 3 .Another example of CEL of a 3-Sat problem with N = 15 and α = 4 is shown in Fig. 3b, when the system is in the Hard-Sat regime, where six local minima are found at E = 1 .More examplar CELs of systems with larger N are found in Fig. S2 of the (SI).In comparison, as shown in Fig. 3c, the CEL of a 3-Sat problem with N = 15 and α = 1 is much simpler in structure without local minimum.This shows that the CELs of small systems show similarity with the expected structure of the corresponding energy landscapes in the "hard" and "easy" regimes in the thermodynamics limit.
In addition, CELs allow us to obtain the statistics of local minima, and the number of local minima n LM is shown as a function of α for the 3-Sat problem in Fig. 2c.As we can see, local minima start to emerge beyond α 2.5 and increase with α .This is again consistent with the phenomenon of increasing algorithmic hardness as α increases.Interestingly, n LM scales with N K−1 , which may imply that the emergence of local minima is related to the number of possible constraints per variable.We further show the distribution of n LM /N K−1 in Fig. 2d, where the distributions become narrower as N increases.We remark that these results are different from most of the previous exhaustic studies on small combinatorial systems which mainly focus on ground states 21 .
We remark that CEL can be readily applied to systems with discrete coupling strength, which include many representative spin models and combinatorial optimization problems.The present approach can also accommodate discrete external fields which commensurate with the magnitude of conpuling strength.Nevertheless, for systems with non-discrete coupling strength and external fields, binning of energy would be a simple way to generalize our method to these systems.

Trapping dynamics
Thanks to the largely reduced number of nodes and the identification of local minima in CELs, they allow us to reveal the complete non-equilibrium dynamics when these glassy systems are trapped in local minima, at an arbitrary temperature for an arbitrarily long time.Here, one can formulate a matrix of transition probabilities T a→b from a cluster a to b, describing the Metropolis-Hasting (Markov Chain Monte Carlo (MCMC)) dynamics of the system following the Boltzmann distribution 22,23 .In this case, T a→b for a = b is given by where E a→b = E b − E a and β is the inverse-temperature; n a corresponds to the size of cluster a, and n a N corresponds to the total number of links connecting its constituent configurations, including those internal links within cluster a.On the other hand, for the system to stay in cluster a, the system can either reject the transition to a configuration outside a or transit to another configuration within a, with a total probability given by T a→a (β) = 1 − b� =a T a→b (β) .We then denote the probabilities to find the system in configurations in individual clusters at time t by a vector � P t = (P 1 , . . ., P C ) , and express where T β is the transition matrix with element T a→b (β).
( www.nature.com/scientificreports/With the matrix T β for specific instances, one can conduct a spectral analysis to compute its eigenvalues, and we denote n (β) to be the n-th largest eigenvalue of T β .We will omit the dependence of n on β in subsequent discussions for clarity.We first diagonalize T β as β by the diagonal matrix � β composed of the eigenvalues of T β , and the matrix Q composed of their corresponding eigenvectors.The power of T β can be computed by β , such that P t in Eq. ( 4) is given by ( 5) www.nature.com/scientificreports/ We remark that � t β can be readily computed by the power of the diagonal element of � β since it is a diagonal matrix.Thus, one can compute the probability of reaching any clusters in the CELs at any time step t.Moreover, by considering the limit � ∞ β = lim t→∞ � t β , one can evaluate the equilibrium probability � P t=∞ of the system taking any configuration clusters after an infinitely long time.
We remark that finding the eigenvalues of T β not only facilitates the computation of P t , but also provides important insights into the structure of the energy landscape and its non-equilibrium dynamics.Since T is a transitional probability matrix, by Perron-Frobenius theorem 24 , the absolute values of its eigenvalues are bounded by 1, and except those eigenvalues equal to 1, all of them vanish as t is large.Therefore, the closeness of an eigenvalue to 1 thus gives us a measure of the metastability the corresponding eigenmode, e.g.trapping in a local minima, and thus how strong the non-equilibrium dynamics is trapped.
As an example, we show in Fig. 4a-c 1 to 10 at different inverse-temperatures β for the spin glass, K-Sat instances with α = 4 and α = 1 shown in Fig. 3a-c respectively.In Fig. 4a, we first note that the eigenvalues of the spin glass instance are in pairs due to the symmetric nature of its energy landscape.More interestingly, as shown in both Fig. 4a and b, n differ more at small β , but the few eigenvalues after 1 start to approach 1 as β increases.The number of eigenvalues approaching 1 is equal to the number of local minima in the corresponding CELs, i.e.
3 and 4 of the spin glass instances correspond to the two symmetric local minima in Fig. 3a and 2 to 7 of the 3-Sat instance with α = 4 correspond to the six local minima with E = 1 in Fig. 3b.The increasing proximity of these eigenvalues to 1 also corresponds to an increasing trapping in local minima when β increases, comparable to the trapping in global minima with 1 = 1 .This also raises a question on whether the systems equilibrate at the ground states at zero temperature (i.e.β → ∞ ), since n → 1 for the local minima and are equivalent to 1 = 1 at the ground states.In comparison, as shown in Fig. 4c, other than the largest eigenvalues 1 , the other eigenvalues of the 3-Sat instance ( α = 1 ) with CEL shown in Fig. 3c do not approach 1 when β increases.Instead, these eigenvalues decrease as β increases.This suggests that these t n 's vanish faster as t grows, implying that the system converges to the global faster at a lower temperature.
Since we can obtain the complete transition matrix for these small systems through CELs, one can compute their complete dynamics at an arbitrary temperature for an arbitrarily long time using Eq. ( 5).Starting with a uniform P 0 , we show the probability P g of the spin glass, and 3-Sat instances with α = 4 and α = 1 being in the ground state after t = 2 17 (blue circles) and 2 37 (red squares) iteration steps as a function of β in Fig. 5.As we can see in Fig. 5a and b, for both spin glass and 3-Sat instances with α = 4 , P g first increases with β as expected, but remarkably decreases as β further increases; the MCMC simulation results are in good agreement with these theoretical predictions by CELs.These results imply that with a random initial condition, the trapping at local minima becomes more significant as temperature decreases below some specific values and it is less likely  The probability P g of finding the ground states in the examplar spin glass, and the 3-Sat problem instances with α = 4 and α = 1 , of which CELs are shown in Fig. 3a-c respectively, as a function of inverse- temperature β obtained by Eq. ( 5) after t = 2 17 and 2 37 iterations, compared with simulation results.Insets: P g as a function of time t.
to locate the ground states within a finite time.Such observations have a strong implication on the cooling procedure in physics-inspired optimization algorithms such as simulated annealing.In comparison, instead of an optimal non-zero temperature, P g for the 3-Sat instance with α = 1 monotonically increases with β , i.e. monotonically decreases with temperature, implying a zero temperature maximizes the probability of locating the ground state as expected.As we can see, time t seems to play an important role as the optimal range of β with high P g widens with t.We further show in the insets of Fig. 5 that P g does increase with t.Nevertheless, for spin glass and the "hard" 3-Sat instance as shown in the insets of Fig. 5a and b respectively, the increase of P g is not smooth and multiple jumps and plateaus are observed, implying that local minima are completing with the global minima for the probability but they cease to trap the system at different time t.This phenomenon can be explained by eigenvalues, where a sufficiently large t makes t n of the local minima sufficiently less than 1, despite n ≈ 1 (see again Fig. 4a and b).This also implies that for glassy systems such as spin glasses and "hard" K-Sat problems, the proximity of n to 1 is related to the ability or stiffness of individual minima in trapping the system, which may depend on their entropy or number of external connections; one may thus estimate the characteristic duration of trapping in individual minima using n .In comparison, for the 3-Sat instance with α = 1 as shown in the inset of Fig. 5c, P g increases smoothly with time t, implying there is trapping in dynamics.
As the eigenvalue n of local minima approaches but is not equal to 1, one may anticipate that ultimately all t n of the local minima would be sufficiently less than 1, and only the global minima left.In otherwords, the MCMC dynamics ultimately lead to the ground states.However, in a practically shorter time-scale, we see that spontaneous ergodicity breaking 25,26 emerges as the time-scale is sufficiently large for the systems to explore the whole configurational space according to the equilibrium distribution.
We remark that MCMC simulations are in good agreement with theoretical predictions, including the drop in P g with β and the abrupt jumps and plateaus of P g at small t in Fig. 5a and b, while the small discrepancies may come from the mean-field nature of the clustered transition probabilities in Eq. ( 5).In addition, since one can easily compute T t β for an arbitrarily large t, e.g. 10 14 in the insets of Fig. 5, by repeatedly powering T t β and its products, one can obtain the long-time dynamics by Eq. ( 5) out of reach by modern computational capability.For the sake of a clear illustration and elaboration, we show the above results for only three instances; in Fig. S3 of the SI, we show that the sample-averaged P g with 200 samples exhibits a similar behavior against β and t.Sample-averaged MCMC simulation results are also in good agreement with theoretical predictions.

Implications on cooling
The eigenvalues of the transition matrix T β from CELs also provide implications on the essence of cooling in identifying low-lying states, e.g.simulated annealing (SA), especially for glassy systems such as spin glasses and "hard K-Sat instances.As we see from Fig. 4a and b, the difference among eigenvalues is large at small β when the system can distinguish 1 which attributes to the global minima from n which attributes to the local minima.As β increases, these largest eigenvalues are getting closer in values.By cooling the system from a high temperature, the system does not start with a random state at the beginning of each cooling stage but instead continuously biases towards the global minima due to difference between its 1 from other n , though this difference is vanishing.This suggests that cooling is more effective in identifying ground states compared with fixedtemperature dynamics.Nevertheless, once the system is trapped in a local minimum, lowering temperature in SA does not help the system escape from the minimum, and only an extremely slow (and potentially infeasible) cooling schedule may help.

Partial coarse-grained energy landscape (PCEL)
The computation of CELs is only feasible for systems with small size N since it requires examining all 2 N variable configurations.Nevertheless, for large systems, we introduce a method to obtain a partial coarse-grained energy landscape (PCEL) for the low-energy configurations.In this case, we sample variable configurations by MCMC simulations at a fixed sampling inverse-temperature β s for T steps, and restart the sampling with random initial conditions for multiple times.We record all the sampled configurations for the construction of PCELs following the same procedures as in CELs.By using an appropriate β s , one can extract specific part of the energy landscape, for instance, a moderately large β s for extracting the low-energy configurations.
We remark that PCELs are only approximations since MCMC simulations with finite time T can only sample a small fraction of all 2 N configurations, though the number of sampled configurations for systems with large N can be significantly larger than those of the small systems we presented before.In addition, it can happen that clusters with the same energy in PCELs indeed belong to a larger cluster since not all the intermediate configurations between the two clusters are sampled.An example of PCEL for a K-Sat problem with N = 50 is shown in Fig. S4a of the SI.Since we are mainly interested in the glassy behaviors contributed by the global and local minima, to further simplify the analyses, we make one more approximation to leave only a single shortest path between minima in PCELs to obtain a simplified transition matrix Tβ ; the simplified version of PCEL in Fig. S4a is shown in Fig. S4b.We found that the results obtained by the simplified Tβ are similar to those without this simplification.
The major advantage in using PCELs to analyze system dynamics is that a single MCMC procedure to extract the PCEL at a single β s can provide us with the dynamics of the system at an arbitrary temperature for an arbitrar- ily long time out of which by simulations.We show the dynamics of a 3-Sat problem with N = 50 in Fig. 6, which is obtained by the simplified PCEL shown in Fig. S4b sampled at a single β s = 5 .The theoretical results agree well with simulations at different β except those at small β when the system explores high-energy configurations while PCELs focus on low-energy configurations.The corresponding sampled-averaged plot is shown in Fig. S5.The same phenomena as in the small systems are observed, namely the drop in P g as temperature decreases, as www.nature.com/scientificreports/well as the jumps in P g as time t increases.These results suggest that the findings based on CELs in small systems are also observed in large systems, which show that CELs and PCELs open up a new platform for us to reveal the long-time non-equilibrium dynamics for glassy systems.

Conclusion
We introduced an approach called Coarse-grained Energy Landscape (CEL) to reveal for the complete energy landscapes of small glassy and non-glassy systems, showing clearly their differences.In terms of methodology, by formulating CELs and analyzing their transition matrix, one can analytically compute the non-equilibrium dynamics of a system at an arbitrary temperature for an arbitrary long time, out of reach by existing theories and numerics, and again revealing clearly the differences between glassy and non-glassy systems.For large systems, we introduce a variant approach to partially reveal the energy landscapes, which allow us to conduct the same analysis as in small systems.In terms of understanding, we a clear and complete physical picture on how glassy systems are trapped in local minima, as evident from the drop in the ground-state probability as temperature decreases as well as their abrupt jumps as time increases.Such phenomena are not observed in non-glassy systems.Simulation results agree well with theoretical predictions.Our work advances methodology by a new tool for analyzing the non-equilibrium dynamics of complex disordered systems, which generate clear, complete and new understandings and insights on their long-time behavior inaccessible by existing approaches.S4, as a function of β , obtained by the transition matrix from PCELs sampled for T = 10 5 steps at β s = 5 with 10 re-starts, and then by Eq. ( 5) after t = 2 17 and 2 37 iterations, compared with simulation results.Insets: P g as a function of time t, where a factor of ln[2 N /n(� s sampled )] has been multiplied to t in the results obtained by PCELs since only part of the energy landscape is extracted, with n( s sampled ) to be the number of sampled configurations.

Figure 1 .
Figure 1.(a) An example of FEL with 2 N = 32 configurations from E = 7 at the top to E = 0 at the bottom, of a 3-Sat problem with N = 5 variables and α = 4 .(b) The corresponding CEL with C = 17 clusters.Global minima and local minima are shown in squares and triangles respectively; node size corresponds to the number of constituent configurations in the clusters; red links correspond to the connections to local minima.

Figure 2 .
Figure 2. (a) The number of clusters in CEL divided by the number of configurations in FEL, i.e.C/2 N , as a function of N. (b) The exponent γ in C ∝ 2 γ N as a function of α in the 3-Sat problem.(c) The number of local minima in CEL denoted as n LM , rescaled with N 2 of the 3-Sat problem as a function of α for different system size N.(d) The rescaled probability distribution P(n LM /N K−1 ) .All the results are obtained by averaging over 10,000, 2000, 1000, 500, 100 realizations for cases with N = 5, 10, 15, 20 and 25 respectively.

Figure 3 .
Figure 3.The low-energy portion of examplar CELs for an instance of (a) spin glass on random regular graph with f + = 0.5 , and 3-Sat problems with (b) α = 4 and (c) α = 1 , all with N = 15 .Global minima and local minima are shown in squares and triangles respectively; node size corresponds to the number of constituent configurations in the clusters; red links correspond to the connections to local minima.

10 Figure 6 .
Figure6.The probability P g for a 3-Sat problem instance with N = 50 of which PCELs are shown in Fig.S4, as a function of β , obtained by the transition matrix from PCELs sampled for T = 10 5 steps at β s = 5 with 10 re-starts, and then by Eq. (5) after t = 2 17 and 2 37 iterations, compared with simulation results.Insets: P g as a function of time t, where a factor of ln[2 N /n(� s sampled )] has been multiplied to t in the results obtained by PCELs since only part of the energy landscape is extracted, with n( s sampled ) to be the number of sampled configurations.